#####################################################################################################
# This script generates two maps of disenfranchisement in antebellum Louisiana, one before          #
# and one after reform.                                                                             #
#                                                                                                   #
# It draws on existing data series on the location of lakes, rivers, canals, and railroad,          #
# as well as the topography of the state.                                                           #
#                                                                                                   #
# Topography data is shown, but not necessary for replication. This can be commented out if you     #
# have difficulties with get_elev_raster.                                                           # 
#                                                                                                   #
#####################################################################################################

### Start by loading packages used across different maps. 
rm(list = ls())

library(shapefiles)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(ggplot2)
library(readxl)
library(dismo)
library(raster)
library(rgeos)
library(elevatr)
library(ggmap) 
library(tibble)

### Provide version list. 
sessionInfo()

### Set the API key to get the Google geocoded location of cities. For users without a key, longitude and latitude are provided below
#register_google('[API key here]')
### Set the directory
dir <-paste("../", sep="")
setwd(dir)

### Create a heat map function for maps. 
plot.heat <- function(states.map,z,title=NULL,breaks=NULL,reverse=FALSE,cex.legend=1,bw=.2,col.vec=NULL,plot.legend=TRUE) {
  if (is.null(breaks)) {
    breaks=
      seq(
        floor(min(states.map@data[,z],na.rm=TRUE)*10)/10
        ,
        ceiling(max(states.map@data[,z],na.rm=TRUE)*10)/10
        ,.1)
  }
  states.map@data$zCat <- cut(states.map@data[,z],breaks,include.lowest=TRUE)
  cutpoints <- levels(states.map@data$zCat)
  if (is.null(col.vec)) col.vec <- color.palette1(length(levels(states.map@data$zCat)))
  if (reverse) {
    cutpointsColors <- rev(col.vec)
  } else {
    cutpointsColors <- col.vec
  }
  levels(states.map@data$zCat) <- cutpointsColors
  plot(states.map,border=gray(.8), lwd=bw,axes = FALSE, las = 1,col=as.character(states.map@data$zCat))
  if (plot.legend) legend("bottomleft", cutpoints, fill = cutpointsColors,bty="n",title=title,cex=cex.legend)
}

### Define color scheme and create index for figures
color.palette1 =  colorRampPalette(c("black","white"),space="rgb", bias=1, interpolate="spline")
cols <-color.palette1(6)
index = c(cols[1], cols[3], cols[4], cols[5], cols[6])

################################################################################
###                        Louisana Census of Voters                         ###
################################################################################

### Set figure details and year range

year <- paste("1841")
#year <- paste("1858")
lgb1		= "10-30%"
lgb2		= "30-50%"
lgb3		= "50-70%"
lgb4		= "70-95%"
lgb5		= "95-100%"

### Set size of New Orleans disenfranchised population depending on year.
if(year=="1841"){
  disf  = 20453
  year2 <-paste("1851")
  label <-paste("Figures/Supplemental/Supplemental-Figure-2-1.pdf", sep="")
} else if(year=="1858"){
  disf  = 26903
  year2 <-paste("1868")
  label <-paste("Figures/Supplemental/Supplemental-Figure-2-2.pdf", sep="")
}

### Create list of inputs for shape/voter files. Shapefiles for Louisiana in 1841 and 1858 (and the state). 
### Voter files for 1841 and 1858.

input1 <- paste("Shapefiles/Louisiana-", year, ".shp", sep="")
input2 <- paste("Data/LA-Census-Of-Voters-", year, ".csv", sep="")
input3 <- paste("Shapefiles/Louisiana.shp", sep="")

### Load rail, rivers, canals, elevation, and the county level data. 
rail <- readOGR("Shapefiles/RR1826-1911Modified0509161", "RR1826-1911Modified050916")
rail <- rail[rail$InOpBy < year,]
rivers <- readOGR("Shapefiles/SteamboatNavigatedRiverContigUSAprojection", layer = "SteamboatNavigatedRiverContigUSAprojection")
rivers <- rivers[rivers$START < year & rivers$END > year2,]
canals <- readOGR("Shapefiles/19thC_Canals_March2017", layer = "19thC_Canals_March2017")
canals <- canals[which(canals$OPENED < year & canals$CLOSED > year),]

counties.map <- readShapeSpatial(input1, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
counties.map <- spTransform(counties.map, CRS(proj4string(rivers)))
state.map <- readShapeSpatial(input3, proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
state.map <- spTransform(state.map, CRS(proj4string(rivers)))
canals <- spTransform(canals, CRS(proj4string(rivers)))
rail <- spTransform(rail, CRS(proj4string(rivers)))
elevation <- get_elev_raster(counties.map, z = 8, expand = 150000)

data.counties <- read.csv(input2)

### The location of the cities requires user to have a Google API key. 
### I comment out the command accessing Google's geocoding, and replace it with the longitute
### and latitude of New Orleans
cities1 <- c("New Orleans, LA")
#cities <- geocode(cities1)
cities <-tibble(lon = -90.0715, lat= 29.9511)
#Input estimated disenfranchised population of cities, from county information and census records.
cities2 <- c(disf)
### Rename New Orleans to parish level (which is where data is from)
cities1 <- c("Orleans Parish")
cities <- SpatialPointsDataFrame(cities, data = cities, proj4string = CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84"))
cities <- spTransform(cities, CRS(proj4string(rivers)))


names(data.counties )[2]<-"FIPS"
names(data.counties )[5]<-"est"

counties.map@data$order <- 1:nrow(counties.map@data)
counties.map@data <- merge(counties.map@data,data.counties, by=c("FIPS"), all.x=TRUE, all.y=FALSE, sort=FALSE)
counties.map@data<-counties.map@data[order(counties.map@data$order),]
names(counties.map@data)[ names(counties.map@data)=="est"]<-"noise"

### Create Supplemental-Figure-2-1 and Supplemental-Figure-2-2 
pdf(label, width=8, height=6, bg="transparent", family="sans" )
par(mfrow=c(1,1),oma=c(0,0,0,0), mar=c(0,0,0,0))

plot(counties.map,density=c(30))
par(new=TRUE)
plot.heat(counties.map,z="noise",breaks=(c(10, 30, 50, 70, 95, 151)), plot.legend=FALSE)
plot(state.map,border=gray(0), add=TRUE)
plot(rail, col = "indianred", add = TRUE)
plot(rivers, col = "dodgerblue", lwd = 1, add = TRUE)
plot(canals, col = "blue", lwd = 1, add = TRUE)

### Create the topographical shading. This is not essential to the replication. 
### If user has difficulty getting elevation data from get_elev_raster above, then
### section can be commented out.
mosaic <- elevation * 50 
slope <- terrain(mosaic, opt="slope", unit='radians')
aspect <- terrain(mosaic, opt="aspect", unit='radians')
hillshade <- hillShade(slope, aspect, angle=45, direction=315)
hillshade2 <- aggregate(hillshade , fact = 5 , method = "bilinear" )
hillshade2 <- focal(hillshade2, w=matrix(1/9, nc=3, nr=3), mean)
slope2 <- aggregate(slope , fact = 5 , method = "bilinear" )
slope2 <- focal(slope2, w=matrix(1/9, nc=3, nr=3), mean)
plot(hillshade2, col=grey.colors(100, start=0, end=1), legend=F, alpha=.1, add=TRUE)
plot(slope2, col=grey.colors(100, start=1, end=0), legend=F, alpha=.1, add=TRUE)

### Plot location of cities / disenfranchised adult free men.
plot(cities, add = T,  pch=21,  bg=alpha(rgb(1,1,1), 0),col="black", lwd=.4, cex=cities2/500)
text(cities, labels=cities1, pos=3, offset=1.5)

legend('topright', fill=index , ncol=1,
       x.intersp = 1, xjust=0, yjust=0,cex=.9, bg="white", box.col = "white",
       legend = c(paste(lgb1, "Adult White Men\nEnfranchised"), lgb2, lgb3, lgb4, lgb5))

dev.off()

